function  [c,f,s] =  bond_pricing_pde( w, lambda_input, v,  dvdw, model_sol, par_bonds, par,  reversal )
% This function generate the necessary results for the PDE solver.
% Refer to https://www.mathworks.com/help/matlab/math/partial-differential-equations.html

if(nargin<8)
   reversal = false;   % note: if reversal = true, then we reverse the definition of lambda as lambda_star - lambda, and solve backward. 
end

% Other values
tau=par_bonds.tau;   pi=par_bonds.pi;    hat_kappa0 = par_bonds.hat_kappa0;   lambda_star=par_bonds.lambda_star;

if(reversal)
    lambda = lambda_star - lambda_input;
else
    lambda = lambda_input;
end

lambda_tau = 1/tau - pi*lambda;
kappap = model_sol.kappap_fitted_matrix(w, lambda);
xK = model_sol.xK_fitted_matrix(w, lambda);
sigmap = model_sol.sigmap_fitted_matrix(w, lambda);
sigmaw = model_sol.sigmaw_fitted_matrix(w, lambda);
kappab = model_sol.kappab_fitted_matrix(w, lambda);
muw =  model_sol.muw_fitted_matrix(w, lambda);
rf = model_sol.rf_fitted_matrix(w, lambda);


% c results
c = par.mu_lambda_fun(lambda) ./ (  w^2.*sigmaw^2 );
f = -1/2*dvdw;
s = 1./( w.^2 * sigmaw.^2 ) .* (  lambda*pi./(1-kappab).*( v - (1-kappap-hat_kappa0) ) + lambda_tau*(v-1) + dvdw.*( xK*(par.sigmaK+sigmap).*w.*sigmaw - w.*muw )  +  rf.*v   );









